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Abstract 



We investigate how stochastic reaction processes are affected by external pertur- 
bations. We describe an extension of the deterministic metaboUc control analysis 
(MCA) to the stochastic regime. We introduce stochastic sensitivities for mean 
and covariance values of reactant concentrations and reaction fluxes and show that 
there exist MCA-like summation theorems among these sensitivities. The summa- 
tion theorems for flux variances and the control distribution of the flux variances is 
shown to depend on the size of the measurement time window (e) within which re- 
action events are counted for measuring a single flux. It is found that the degree of 
the e-dependency can become significant for processes involving multi-time-scale 
dynamics and is estimated by introducing a new measure of time scale separation. 
This e-dependency is shown to be closely related to the power-law scaling observed 
in flux fluctuations in various complex networks. We also propose a systematic 
way to control fluctuations of reactant concentrations while minimizing changes 
in mean concentration levels by applying the stochastic sensitivities. Key words: 
metabolic control analysis; sensitivity; stochastic process; noise propagation 
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Introduction 

Metabolic control analysis (MCA) ([l], Q B and the closely related biochemical 
systems theory have greatly enhanced our ability to understand the dynamics of 
cellular networks. However, these approaches are based on a deterministic picture 
of cellular processes and in recent years it has become clear that many networks, 
such asgene regulatory networks, operate with a significant degree of stochastic- 
ity fl la B, HJa, fiob . In these situations a deterministic formalism is inadequate 
dllL I12I. I13I . Il4h . In this paper we begin the process of developing a new analy- 
sis method of control on stochastic dynamics by extending MCA to the stochastic 
regime. We call the extension stochastic control analysis (SCA). 

MCA is an analysis of sensitivities which quantifies how much system vari- 
ables change in response of the perturbations in system parameters. To extend 
the MCA to the stochastic regime, we need to introduce sensitivity measures for 
stochastic system variables. There have been a wide variety of efforts in recent 
years to introduce and investigate sensitivity measures for stochastic systems re- 



lated to mean levels of concentrations and fluxes in stochastic reaction systems (111 
|20|) . More pertinent to this paper is the work by Andrea Rocco 
who investigated the MCA summation and connectivity theorems related to the 
most-probable concentration values and their con^esponding reaction rates (i21i) . 
However, the sensitivities for noise characteristics (variance and covaraince) were 
not investigated and the summation theorems related to the noise properties were 
not discussed. Thus, a systematic MCA-like approach for controlling noise has not 
been made. 

In this paper we will focus on the control coefficients for variances 

and covariances of concentrations and fluxes. The control coefficients quantify the 
global responses due to (static) perturbations in the system parameters. We also 
introduce sensitivities for the mean levels of concentrations and fluxes, which are 
closely related to the MCA control coefficients. We obtain MCA-like summation 
theorems for the stochastic sensitivity measures. In a similar way to the determin- 
istic MCA theorems, the SCA theorems imply that control is distributed over a 
reaction system. 

The summation theorems for flux variances show very interesting properties: 
Flux is measured by counting the number of reaction events within a given time 
window e. From this we show that the sum value can be highly dependent on 
the measurement time window (e). This in turn implies that the control of flux 
variances can be sensitive to the value of e. We present a case where a control dis- 
tribution over a reaction system changes significantly with the value of e and thus 
the system can show an increase or decrease in the flux variances depending on the 
value of e. The degree of such e-dependency of the flux variance control is closely 
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related to how both the time scales of the fast and slow fluctuations are separated. 
Such a separation can be quantified by introducing a new time scale separation 
measure which can be estimated from the temporal sequences of reaction events. 

The summation theorems for flux variances also show a close connection to 
the scaling relationship between flux variances and their mean values recently 
observed in various complex networks: the Internet, microprocessor logic net- 
works, the World Wide Web, highway systems, river networks, and stock markets 
(I22I. I23I |24|. I25I. l26h. In these systems, fluxes were defined as the number of packets 
processed in network routers in the Internet, activity of connections between logic 
gates in the microprocessors, the number of visits on sites in the World Wide Web, 
the number of cars of traffic at different locations in the highway systems, stream 
flows in the river networks, and the traded values of stocks in the stock market. It 
was investigated how the standard deviation (cr) of the flux is related to the mean 
value of the flux ((/)): cr ~ (/)"■ de Menezes and Barabasi claimed that the Inter- 
net and microprocessor logic networks belong to a universality class characterized 
by an exponent value of a = 0.5, and the World Wide Web, highway systems, and 
river networks to that of a = 1 (22). However, stock markets such as NYSE and 
NASDAQ show non-universal values of a (|24l. l25n. Meloni et al. (\26i ) proposed a 
model of random diffusion to show how the value of the exponent can crossover 
from 0.5 to 1 depending on the number of links connected to a node, the strength 
of external noise and the time measurement window size e. In this paper, we show 
a connection between the summation theorems for flux variances and the scaling 
crossover phenomena. We briefly discuss that the exponent crossover can take two 
different forms depending on the time window size e relative to the coiTclation time 
of the external noise. 

As an application of the stochastic control coefficients, we provide a systematic 
non-local method for the control of noise levels in concentrations, while minimiz- 
ing changes in mean concentration values. Such orthogonal control is performed 
first by estimating control coefficients for mean values and coefficients of varia- 
tion (CVs) of concentrations under the linear noise approximation ilm . From the 
estimate of the control coefficients, we find the direction of the parameter pertur- 
bations leading to a sensitive response of change in the CVs of the concentrations 
while minimizing the change in their mean values and this enables us to identify 
which sets of parameters need to be controlled by how much in a relative sense. We 
apply this orthogonal control in a negative feedback system under external noise 
and successfully reduce the concentration noise level, with a negligible change in 
the concentration mean value. 
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Model Systems and Definitions of Control Coefficients 

The model system we will consider is a chemical reaction system described by the 



chemical master equation (1281. 1291). i.e., we assume the system is spatially homo- 
geneous (uniform concentrations throughout the time evolution of the system). We 
assume that it can be described by L kinds of reaction rates for M reactants. The 
system is composed of the external and internal processes. The external process is 
modeled by allowing one of the species (denoted by either Se or Si) to fluctuate 
slowly and independently, compared to the rest. Se is considered a source of ex- 
ternal noise. The internal system, composed of all other species, is affected by the 
external noise and also by internal noise generated from the internal reactions. 

To estimate how a system responds under parameter perturbations at the sta- 
tionary state, we introduce sensitivity measures called control coefficients. The 
system variables (y) of interest can be either mean values or coefficients of varia- 
tion/covariation (CV/CCV) of concentrations and reaction fluxes. CV is variance 
divided by the mean square and CCV is the covariance (between two variables) 
divided by the product of their mean values. We define the control coefficients for 
these variables as 

(jy ^ ^ dlogy 
^ y dp d log p ' 

which indicates the relative change in y due to a given relative change in a param- 
eter p. The change in y is from one stationary state to another corresponding to 
before and after the perturbation, respectively. We note that control coefficients for 
different system variables - most-probable concentrations (not mean concentra- 
tions) - have been investigated in the framework of MCA, but sensitivities related 



to fluctuation properties have not (1211) . The parameter p will be called here a con- 
trol parameter, which is not affected by the system's reactions. We restrict the set 
of the control parameters (p = (pi, • • • ,pl)) to be the proportionality constants of 
reaction rates. E.g., for a reaction rate v = with s concentration and Km 

a Michaelis constant, p is a control parameter but Km is not. The total enzyme 
concentration that catalyses a reaction is one such parameter. 



SCA: Summation Theorems for Control Coefficients 

We have found that there exist MCA-like summation theorems among the proposed 
stochastic sensitivities, which are valid under any strength of noise and finite per- 
turbations of parameters p. The existence of these theorems is rooted in the fact 
that the stochastic measures satisfy certain scaling properties under a specific kind 
of scale change in time and control parameters. 
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Summation theorems for concentrations 

We note that all reaction propensity functions Vi are proportional to control pa- 
rameter Pi'. Vi{s,ap) = avi{s,p). Let us scale all control parameters by a fixed 
proportion a. The simultaneous change in all propensity functions can be inter- 
preted as a change in the time scale in the amount of 1/a because the propensity 
functions are inversely proportional to time. Mean levels, CVs and CCVs of con- 
centrations are time independent variables at stationary states. This means that 
these quantities remain the same under the parameter change Jioh . We can sum- 
marize these arguments with the following equation (refer to Table[T]for notation). 
The change in a concentration mean level is expressed as: 

Pi 

for all j = 0, • • • , M. Since d{sj) = 0, we derive 

L 

1=1 

for all species j. The same argument can be applied for the concentration CVs and 
CCVs. 

L 

1=1 

for all species j and k. 



Summation theorems for fluxes 

To derive the summation theorems for mean fluxes, we consider again the param- 
eter scale change. Under this change, mean propensity functions will scale by a. 
Since the mean propensity function {vi) is equal to the mean fluxes {Ji), the mean 
fluxes will also scale by a (|3d) . Since the scale change in the mean flux can be 
expressed as: 



{Ji) fr( Pi P ^ 

we obtain summation theorems for mean flux control coefficients: 

L 



^Cj,^'^ = l. (3) 
1=1 
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We will also derive summation theorems for flux CVs and CCVs. However 
before we derive them, it is important to clarify the difference between a propensity 
function, a reaction rate, and a reaction flux. All of them are stochastic variables. 
The reaction flux J is measured by counting the number of reaction events within 
a time window e: 



Ji 



Number of events of (i-th kind) reaction 
occurred during e 



The propensity function is a mathematical function previously denoted by v. The 
mean values of both v and J are equal. The fluctuation strengths of each can 
however be different, because the variances of J are dependent on e (as will be 
discussed later), while those of v are not. We express the CV/CCV of J {V'^) 
as a function of e: V'^{e,p). We use the term, reaction rate, as either the flux or 
propensity function, depending on context. 

Now we will derive the summation theorems for flux CVs. The first thing to 
note is that flux CVs are unitless in time. The flux CVs obtained by scaling all 
parameters by a is the same as those obtained by scaling the time by 1/a: 

V-^ {e,ap) = W {ae,p). 

This can be expressed in terms of control coefficients as follows: 

Y.Cvt = ^^, (4) 



=1 



for all reactions j, k. This equation means that the sum value is equal to the slope 
of a log-log plot of flux CV and CCV vs. e. Since the flux CV and CCV depend on 
e, the sum value can also depend on e. 



Summation theorems for flux CVs vs. multi-time-scale dynamics 

In this section, we investigate how the sum value of Eq.|4] changes with e. We have 
found an interesting fact that the sum value can vary significantly with the change 
in e when the system shows wide distributions of reaction time scales. 

We briefly discuss the mechanisms for leading to the sum value change by 
considering a simple reaction system: a two-step reaction cascade as shown in 
Fig- [IK- 5*1 is created with a rate vi and degrades with a rate V2- Si enhances the 
conversion of X2 to ^2- We assume that the creation and degradation processes 
of Si are much slower than those of ^2- Si is the source of external noise. The 
reaction process involving 5*2 is considered an internal system. The time evolution 
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trajectory of 5*2 shows a mixture of two different kinds of noise (slow and fast) as 
shown in Fig. [1^ Jiil. 23, 2^. In the time resolution of e 0.1, external noise 



is negligible while the internal noise (caused by reaction events of ^3 and ^4) is 
dominant. As e increases, the external noise becomes more dominant while the 
internal noise becomes averaged out. The creation flux of S2 shows this tendency 
clearly as shown in Fig.[T]Zl. 

We have plotted all internal and external flux CVs and CCVs vs. e (Fig. [TjD). 
First, we discuss the flux CVs: K/. For the fluxes con^esponding to the fast reac- 
tions (i = 3, 4), a plateau region appears (slope 0, i.e., the sum value ~ 0) and for 
the fluxes corresponding to the slow reactions {i = I, 2), they don't. The plateau 
region appears due to the fact that the internal noise becomes sufficiently averaged 
out at the time scale of e ~ l/pa = 1 and the external noise becomes dominant 
for all values of e > l/ps- J3 can be approximated to be ^3 for e 10 (Fig.[T]C): 
J3 ~ 773 = ^351. The flux CV can be expressed as V33 ~ Vf^ = = 0.1. This 
is the value of the flux CV at the plateau region. This is because the external noise 
has the correlation time (I/P2) and the approximate equality J3 ~ V3 persists until 
the external noise is coiTclated in time, up to e ~ l/p2 = 100. Therefore, the slope 
of the log-log plot of V'^ vs. e becomes close to zero (Fig. |2l), which means that 
the sum value of the flux CV is also close to zero. 

For e ^ r(= I/P2), Si does not fluctuate compared to 82- S2 can be consid- 
ered to be created from a constant source. The probability P(n; e) of having the 
number n of events of reaction during time e satisfies a Poisson distribution: 

P(n;e) = e--i^. 

n! 

The flux CV becomes inversely proportional to e (Fig.|2l): 



J {Ji)-{J3)' in') -in)' 1 



33 / 7„\2 



(J3)2 (n)2 (n) e(J3) 

Thus, the sum value of the flux CV control coefficients is -1. 

For e » r, the external noise becomes uncorrelated in time at this time scale. 
Thus, the flux J3 measured by using this e value will be uncorrelated (statistically 
independent) in time. We denote the minimum of such a value of e by eind- For the 
value of e » eind, the flux estimate J'^ can be considered an average of independent 
samples of J'^»"d with a sample size el eind- Therefore, V'^^ = y>^^™'' 1 — qc -. 
(From Fig. |2l e^nd is ~ 10^.) This explains intuitively why the flux CV scales as 
1/e for large e values (Fig. |2l). Therefore, the sum value of the flux CV control 
coefficients is - 1 . 
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For each different pair of fluxes, the asymptotic form of its coefficient of co- 
variance for e ^ r, is different: either a plateau or a straight line proportional to e. 
A detailed discussion on this is provided in the Supporting Material. 

The arguments presented above can be generalized for a typical reaction sys- 
tems showing flux fluctuations with two different time scale dynamics (refer to the 
Supporting Material). A plateau region (for intermediate e) and two regions of -1 
slope (for very small e and large e) can appear typically for CVs of such fluctua- 
tions. 

Estimation of time scale separation 

As presented previously, the plateau region in Fig. |2] appears due to the time scale 
separation between fast and slow system dynamics. If the separation is not wide 
enough, the plateau region can be tilted. In this case, the sum value of the flux CV 
control coefficients will deviate from zero in the region of the plateau. To identify 
such deviations, we propose a time-scale separation measure. 

The separation measure ($) quantifies the vertical distance between the two 
asymptotic linear lines for the log-log plot of flux CV vs. e coiTcsponding to e — > 
and oo as shown in Fig. [2] The larger the measure the wider the plateau region 
and the smaller its slope, i.e., the sum value of flux CV control coefficients becomes 
closer to zero. Consider a reaction step with its propensity function given by v{se), 
explicitly showing a dependency on the external species Se acting as a source of 
slow noise. If events in the above reaction t^(se) do not cause lai^ge fluctuations in 
substrate concentrations, we can propose the separation measure to be: 



where A is the area underneath the auto-coiTclation function of the external noise: 



The area A can be further simplified to be the variance of Sg multiplied by its cor- 
relation time (r) as a first level of approximation. The derivation of Eq.|5]is given 
in the Supporting Material. The measure increases with an increase in either the 
sensitivity of the propensity function (dv/ds), the absolute strength of the exter- 
nal noise, or the correlation time of the external noise. The measure, however, 
decreases with an increase in mean flux (f (se)) (Fig-EJl. This is counter-intuitive 
because: the larger the flux, the faster the concentration fluctuations and the wider 
the time scale separation. However, the increase in the mean flux, depending on 
which parameters to control, can lead to an increase in the time scale separation via 




(5) 



A = lim 
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the change in the sensitivity. E.g., If global proportionality constants are increased 
by x%, both the sensitivity and mean flux increase by x%. Thus, as a net effect, 
the measure can increase for this choice of control. 

Power-law scaling in flux fluctuations 

We will now show briefly how the slope change is related to power law sc aling that 



is observed in flow fluctuations in other complex networks (1221. 123l. I24l |25|, 1261) . 
In the scaUng studies, it was investigated how the flux CV is related to mean flux 
(actually, rather than flux, but the number of events occurred within e, was investi- 
gated). As shown in the Supporting Material, depending on the value of e relative 
to the correlation time of the external noise, the scaling crossover takes different 
forms (Eqs. S2 and S3). We propose here that the scaling crossover that appears 
in other complex networks can also depend on the interplay between the external 
noise correlation time and e. We note that only in the case of a linear propensity 
function v(s p) = ase for e <C r, we could regenerate the crossover function given 



in Eq.(7) in (l26h (here the relative noise strength is given by Vaiiance(se)/(se)^> 
and Eq. S4 is used.) We have therefore shown a connection between power-law 
scaling and flux fluctuations in reaction networks. 



SCA: Parametric Control of Noise Level 



In the literature on deterministic control theory (1311) and MCA (1321. 13 3L 1341 1351) 
some authors have considered the orthogonal control of system variables such as 
flux and species concenti'ations. Here we consider the orthogonal control of mean 
concentration levels and concentration CVs, in order to control noise independently 
of the mean concentration levels. Such control needs to satisfy the following re- 
quirements. First, the concentration CV decreases as the concentration mean in- 
creases, and thus the control of mean and CV can be strongly anti-correlated. In 
this case one needs to find systematically which parameters to be perturbed by 
how much for the orthogonal control. Second, the concentration CV is dependent 



on noise propagation (1421 . 1431) . implying that a set of multiple parameters may need 
to be controlled simultaneously to achieve a sensitive change in CV. Taking into ac- 
count these requirements, we present a systematic non-local method for orthogonal 
control using the control coefficients. 

We introduce a control vector Cp = {Cp-^, Cpj, • • • , Cpj^) defined in an L- 
dimensional control parameter space. When parameters p are perturbed in the 
direction of Cp, a system variable y (concentration mean or CV) shows a sensitive 
response of increase. When p are perturbed in one of the perpendicular directions 
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of Cp, the system variable y does not change. 

We aim to find parameter perturbations (A) that lead to a decrease in the con- 
centration CV without changing the concentration mean. For the mean (s) not to 
be changed, the parameters must be perturbed in the perpendicular directions of 
Cp ^ (Fig. 131). It needs to be determined which one of the perpendicular directions 
leads to a sensitive response of a decrease in the concentration CV. This determi- 
nation can be done by projecting Cp° onto the parameter space perpendicular to 

(s) 

Cp and multiplying by —1: 



|Cp I 



(6) 



where the factor of — 1 appears since should decrease. 

The efficiency of this orthogonal control can be estimated by how much per- 
centage ratio of the control vector for CV is projected onto the perpendicular- space: 
I sin where 6 is the angle between the two control vectors. If 9 is close to —180°, 
the two controls are anti-correlated and the efficiency is 0. If 9 is close to 90°, 
the two controls are already orthogonal and the efficiency is ~ 1. 

We provide an example of orthogonal control to reduce the concentration CV 
by investigating a linear chain reaction system (Fig.|5K)- This system is under neg- 
ative feedback control and receives external noise via 5*1 . We can predict the distri- 
bution of control based on the linear noise approximation (Fig.|5^). We have esti- 
mated control vectors for the concentration mean and CV, C^*^ and C^'" (Fig.|5^, 
crosses), by using the Lyapunov equation (refer to the Supporting Material) (also 



known as the fluctuation dissipation relationship (l4lU42l) ): 



Ja + <T^J^ + D = 0, (7) 

with J the Jacobian matrix, a concentration covariance matrix, and D diffusion 
matrix. We estimated 9 and A for the original parameter values. We perturbed the 
parameters along A and estimated the new 9 and A for the perturbed ones. After 
two more iterative perturbations, we could reduce the noise level by 25% without 
changing the mean level (Fig.[5tl!). 



SCA: Flux Fluctuation Control 

In this section, we discuss a way to reduce the flux CV. Consider a scenario where 
a metabolic engineer aims to reduce the fluctuations in the production rate of an 
end product. To this aim, her/his first guess is that reducing the concentration fluc- 
tuations will lead to a reduction in the rate fluctuations. The engineer introduces a 
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negative feedback to reduce the concentration fluctuations. The question we might 
ask is whether this operation guarantees that the rate fluctuations are reduced? 

Let us consider the previous example of the linear-topology reaction system 
with negative feedback (Fig.[5]\). We aim to reduce the fluctuations of Jg, by con- 
trolling p6- Based on Fig.|5^, decreasing pQ causes a reduction in the concentration 
CV of Si. We can decide to reduce pQ to decrease the flux fluctuations. To confirm 
that, we have estimated its flux control coefficients based on stochastic simula- 
tions. We found that the sign of the control coefficients Cp^^ is however negative 
for e < T/ (t/ ~ 1: feedback time scale) and positive for e > tj. Reduction of the 
concentration CV causes an increase in the flux CV for e < tj, while it does not 
for e > Tf. This means that controlling p^ can have an opposite effect depending 
on e. This is due to the fact that flux fluctuations become dominated by different 
sources of noise depending on e. Therefore, in this case, we need to choose the 
appropriate value of e depending on the rate fluctuations caused by which source 
of noise to be reduced. 

A question that comes next is: why does the control distribution of flux CVs 
change with the value of e? Consider again the negative feedback system. There are 
three different time scales, related to the internal turn-over reactions (tj), feedback 
controls (tj), and external noise (r = l/p2)- For our choice of parameter values, 
Ti < Tf < T. Depending on where the value of e resides, the flux CV control takes 
different distributions as shown in Fig. |6] 

Foi' Tf ^ ^ ^ the control distribution for downstream flux CVs is quite 
similar to that for the CV of ^4 (Figure[6^ (e = 40) is compared with Fig.[5^). This 
is due to the strong negative feedback and the slow fluctuation components of 54. 
The reaction flux J3 has been confirmed to be approximately equal to f 3 (-Si, 51) 
with ^4 the slow component of fluctuations of ^4 (graph not shown). 

For e > T, the external noise becomes averaged out. Thus, the downstream flux 
CV becomes less sensitive to the external noise correlation time t(= 1/^2)- That 
is why the control by p2 becomes weaker (see Fig. |6^(e = 1000)). This change 
leads to the change in the sum value of control coefficients for downstream-flux 
CV, from ~ to -1. 

For e <C T/-, the internal noise becomes dominant. All the downstream-flux 
CVs asymptotically follow l/e( J) with (J) the downstream-flux mean. Therefore, 
the control vector for the downstream-flux CV is completely anti-parallel with that 
of (J) (Fig. |6^, e = 0.01), implying that orthogonal control is impossible. 

For e ~ Tf, there is a hump in the plot of flux CV vs. e. This is due to the 
negative feedback, where fluctuations of ^4 can be fed back at the same time scale 
without losing its control strength (see Fig. |6K)- 
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Conclusion 

In this paper we describe extensions of metabolic control analysis into the stochas- 
tic regime for general biochemical reaction networks. We have shown that there 
exist MCA-llke summation theorems for stochastic sensitivity measures for mean 
values and coefficients of variation/covariation (CV/CCV) for concentrations and 
reaction fluxes. The summation theorems for the reaction fluxes have shown that 
the sum values of control coefficients for flux CVs/CCVs depend on the size of the 
measurement time window (e). Such dependency becomes stronger as the reaction 
systems shows multi-time-scale dynamics, i.e. the time-scale separation between 
slow and fast modes becomes larger. We have provided a measure to quantify such 
separation. 

In terms of the stochastic sensitivity measures, we have provided a non-local 
systematic way to control mean values and CVs of concentrations orthogonally. 
We hope this method will be useful for controlling noise levels in various reac- 
tion networks such as gene regulatory networks, metaboUc reaction networks, and 
protein-protein interaction networks. 

Finally, we have shown that the control distribution of flux fluctuations can be 
significantly different depending on e which reflects the dynamics at different time 
scales that emerge under varying values of e. Depending on which noise source the 
flux fluctuations to control is caused by, the appropriate window size needs to be 
chosen. 

This work was supported by a National Science Foundation (NSF) Grant in 
Theoretical Biology 0827592. Preliminary studies were supported by funds from 
NSF FIBR 0527023. The authors acknowledge useful discussions with Hong Qian. 
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{f){x) 


Ensemble average of / (over x) 


at a stationary state 


p 


Control Parameter 


S 


Concentration 


V 


Reaction propensity function 


J 


Reaction flux 




Coefficient of co- variation (CCV) between i and j 


Vn 


Coefficient of variation (CV) of j 


Table 1: Notation 
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Figure Legends 

Figure [B 

Two step cascade reaction system: 5i down-regulates the reaction creating 52 (A). 
The reaction rates involving Si are. set 100 times slower than those involving 82- 
Si applies an external noise onto the (internal) system of 82- Time evolution of 
Si and S2 is shown (B). The region of t = [100, 120] is expanded (B,bottom). 
The time evolution profile of ^2 follows the external noise with rapidly fluctuating 
internal noise (B,top). In the time scale of the order of 1, ^2 does not fluctuate but 
Si fluctuates significantly, i.e., the internal noise becomes dominant (B,bottom). 
J3 is measured with three different time window sizes, e = 0.0625, 8, 1024 (C). 
J3 matches with for e ~ 8, because the internal noise is averaged out, i.e., the 
external noise is dominant in this time scale (C,middle). Flux variance of J3 de- 
creases with the time window size e (C,D). V33 shows a plateau, while Vy^ does not 
(D) {V22 overlaps with VA, and V44 with V33 [not shown in graph]). The stochas- 
tic simulation algorithm (1391) was used. Parameters: {Xi, X2,Pi,P2,P3,Pa) = 
(1,1,0.1,0.01,1,1). 

Figure 121 

Two step cascade reaction system (Fig. [TK): The estimate of V33 from the simula- 
tions is compared with the exact analytic result (Eq. S4) and its asymptotic forms 
corresponding toe^r(=l/p2) and e » r. The two asymptotic lines have slope 
-1 and their vertical separation (at the log-log scale) is denoted by <I>, the time-scale 
sepai^ation measure (Eq.|5]). 

Figure m 

Time-scale sepai^ation measure <^ (Eq. [5]l is verified with numerical simulations. 
External noise Sg, generated by Xq ^^^°> Se ^'^^^> 0, is applied onto a reac- 
tion: v{se) = P3 + k'^+s^ ■ "^^^ reaction flux of v{se) is numeri- 
cally estimated by using Eq. SI (solid line. Original). We have reduced <1> by 
perturbing one or more factors affecting $ for each case (other solid lines). We 
normalized V'^ such that its normalized value for e = 0.1 equals 1 for ease of 
comparison. <I> given by Eq. |5] is shown to predict the separation well (dotted 
lines: log (Normalized V'^) = $ — 1 — log^Q(e)). Parameters used: Xq = 1 for 
all the cases, {pi,P2,P3,P4, K.m,n) = (0.2,0.01,0,100,400,2) for "Original", 
(0.2, 0.01, 0, 100, 20, 1) for "dv/ds", (0.4, 0.01, 0, 100, 400, 2) for "A", (0.2, 0.01, 100, 100, 400, 2) 
for "(v)", and (0.4, 0.02, 0, 100, 400, 2) for "r". 
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Figure 131 

Control vectors Cp = {Cp^,Cp^, ■ ■ ■ ,Cp^) for a mean concentration, (s), and 
its CV, y*, are shown in an L-dimensional parameter space as respective green 
and purple arrows, oriented in different directions {9: the angle in between). The 
control vector Cp" is projected onto a space perpendicular to Cp ^ (turquoise blue 
plane). The projected vector of Cp" is denoted by —A. When pai^ameters pertur- 
bations are directed along A, the CV shows a sensitive response of decrease while 
mean concentration remains the same. 

Figure H 

Orthogonal control of concentration CV for a linear topology reaction system with 
a negative feedback and under external noise (A). Distributions of (scaled) con- 
trol coefficients for mean values and CVs of all species ai^e estimated by perturbing 
each parameter by 5% (B) (simulation: hash and solid bars, linear noise approxima- 
tion [shown only for S4]: crosses). The original parameter set is (Xi, X2, Km) = 
(l,l,104),Po = (Pi,P2,P3,P4,P5,P6) = (0.1,0.01,2 X 104,1,1,1). We aim to 
reduce the noise level of ^4 without changing its mean. Control vectors for {S4) 
and are estimated from the linear noise approximation. Then, we estimated 9 ~ 
164° and A (0.0, 0.0, -0.1, 0.1, 0.1, -0.1). We perturbed Pq in the direction of 
A by 40%, i.e., 6pQ = (0, 0, -8 x 10^ , 0.4, 0.4, -0.4) and the new values of param- 
eters ai-e set to po + ^Po = (0.1, 0.01, 1.2 x 10^,1.4,1.4,0.6). We repeated this 
procedure twice more. The probability distribution functions of 54 are shown for 
the series of the perturbations (C). The final p is (0.1, 0.01, 5000, 2.6, 2.6, 0.24). 
Efficiencies (| sin(0) |) of orthogonal controls are shown for different values of a pa- 
rameter pq for the other parameters fixed (D). Orthogonal control is most efficient 
around 0.2 < pe ^ 0.6. 

Figure m 

Flux control distributions for the linear-topology reaction system with negative 
feedback (Fig. (IK)- Flux CV of internal fluxes (A,left) shows humps in the time 
scale of the feedback (e 1); correlation functions between S4, and all inter- 
nal species are shown on the right. G^y{At) = \imtg^oo{Sx{to + At)Sy{to)) — 
{Sx{to)){Sy{to)). The pai^ameter set pg is used (Fig. [5]caption.) Control distribu- 
tions (control vectors) for flux CVs are significantly dependent on e. For e = 0.01, 
the control vector for each mean flux is anti-parallel with that for each correspond- 
ing flux CV. For e = 40, control distribution becomes similar to the control distri- 
butions of CV of 54 (Fig. 12^.) 
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1 Time scale separation measure 

We consider a reaction step with its propensity function given by v{se), showing a 
dependence on the external noise. If events in v{se) do not cause large fluctuations 
of its substrate concentrations, the counting process of the reaction events can be 
described by a doubly-stochastic Poisson process jicj). The probability P{n; e) of 
having the number n of events of reaction v within the time window e is given by 



where P{{se{t)}) is the probability of having a trajectory of {se{t)} for the region 
t e [0, e], and N{e) = dt v{se{t)). The CV of its flux is given by 



Systems 
(Supporting Material) 



Kyung Hyuk Kim and Herbert M. Sauro 




(8) 



where 



Variance (iV(e)) 




with Sv{Se) = v{Se) - {v{Se))s,- 

For e ^ T, Variance(A^(e)) can be simplified as 




Thus, we obtain 




(9) 



For e ^ r, Variance(A^(e)) is simplified as 
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where Gv{t') denotes an autocorrelation function defined by {Sv{tQ + t')6v{to))s^- 
This can be further simplified by 2eA', with A' = dt'Ggy(^g^^ {t'). A' is the area 
underneath the autocorrelation function. Thus, we obtain 



J 1 / 2A'\ 



for e ^ r 



(10) 



If the fluctuations in Se is mostly confined to the linear region of v{se), then 
6v{se) — a5se with a = '^^^j^^ |g=(se)- Thus, we obtain two different asymptotic 
forms of the flux CV for e ^ r and e ^ r, respectively: 



1 



e{J) 



1 + e 



dv 

dSe 



2 Variance fsp 



and 



e{J) 



1 + 2 



{J) 
dv\'i A 



for e < T, 



dsj (J) 



for e » r, 



(11) 



(12) 



where A = dt'Gss^{t') is the area underneath of the autocorrelation function 
of the concentration of the external species, Sg. 

The time scale separation measure is defined as the vertical distance between 
the two asymptotic linear- lines for the log-log plot of flux CV vs. e corresponding 
to e ^ and oo. The measure is obtained from Eq.dTTI) and ([T2l ): 



^ = log 



^/'dv\'^ A 

(J) 



For the two-step reaction process as shown in Fig. 1 A, we can obtain the fol- 
lowing exact result from Eq.® without any approximation: 



1 



33 



+ 



2 x-l + e- 



e{J3) {Si) 



(13) 



where x = P2^- The second term converges to 1/2 for x ^ 1 and vanishes as 
1/x for X ^ 1- While control parameters are fixed, we vary the value of e (see 
Fig.2). V33 ~ l/e(J3) for e ^ l/p2- As e increases, the flux variance reaches a 
plateau region following l/e(J3) + As e » ^/P2, it follows (1/(J3) + 

2/(S'i)p2)/e- The time scale separation measure for this system is expressed as 



$ = loe 



1 + 



2(J3) 



The measure increases with the internal flux (J3) and the time scale separation 
becomes larger as the internal dynamics becomes faster. 
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Figure 7: Flux covariances of different pairs of reactions in the two step cascade re- 
action system Fig. 1 A. Parameters: {Xi, X2,pi,P2-,P3,P4:) = (1,1,0.1,0.01,1,1). 



2 Coefficients of covariance of fluxes vs. e 

In this section, we will investigate how the sum value of Eq.[4] changes with e for 
the coefficients of covariation(CCV) between two different fluxes by investigating 
the slope of the log-log plot of flux CCV vs. e. For ease of presentation, we 
will consider covariances of fluxes rather than the coefficients of covariation. A 
covariance between two different fluxes is defined as 

4 = {{Jr - {Ji)){J, - (Jj))) = {J^Jl) - {Jr){Jj). 

where the ensemble average (.) is performed over the stationary states obtained by 
independent runs of stochastic simulations. 

Consider a two-step cascade reaction system as shown in Fig.lA. First, we will 
investigate how the flux covariance behaves in the limit of e — > 0. Flux covariances 
show different asymptotic behaviors in the limit of e — > depending on the differ- 
ent pairs of fluxes (see Figj7]l. We will explain the mechanisms that generate the 
different behaviors. 

First, we investigate the flux covariance between Ji and J2. If we assume 
that Ji and J2 become independent in the limit of e ^ 0, the covariance cr/2 
vanishes. This, however, is not what we observed by simulation. This indicates 
that there is a con^elation between them. The coiTclation is due to the fact that one 
reaction of vi will increase Si by one, resulting in the increase of V2 and affecting 
the proabability that the reaction f 2 will occur. We take into account this causal 
correlation to estiamte the flux covariance. For a sufficiently small value of e, the 
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Figure 8: Flux co variance of two reactions vi and V2 in the two step cascade reac- 
tion system Fig. 1 A. Parameters: (Xi, X2,pi,P2,P3,P4) = (1,1,0.1,0.01,1,1). 

dominant contributions to the flux covariance come from two cases: first, reactions 
of vi and V2 occur once for each within the time interval e, with the reaction vi first 
and then the reaction V2, and second, each reaction occurs in the opposite order. 
The contribution of the first case to the estimation of (Ji J2) is, for sufficiently 
small e, 

The contribution of the second case is 

{viV2{Si)) 



Thus, we obtain the covariance: 



J {viV2{Sl + I) + ViV2{Sl)) 
0"12 - 7, 



{Vi){v2). 



Since vi is constant {vi = pi) and V2 = P2S1, we obtain 



1 



1 



0"12 - 2V1P2 = -P1X1P2. 

We have verified this result with the simulation data as shown in Fig. [8] 

The flux covariance between Ji and J3 in the limit of e ^ can be also 
estimated in the same way as above: 



4s 



1 



1 



-VlP3 



-piXips. 
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Figure 9: Flux co variance of two reactions vi and in the two step cascade reac- 
tion system Fig. 1 A. Parameters: (Xi, X2,pi,P2,P3,P4) = (1,1,0.1,0.01,1,1). 



The covariance is estimated at 0.05 (see Fig©. 

£7/4 converges to linearly with e as e ^ 0. This is because an event of 
reaction vi does not make any change in the number of 52. The only way to make 
a coiTclation between Ji and J4 is thi^ough an event of reaction v-^. By taking into 
account such indirect effects, the contribution to (Ji J4) becomes 

^ dt dt' £ dfviVsiSi + l)vi{S2 + l)P(5i, 52) 



pmpMSi + 1)(S2 + l)>e. 



Since the non-zero effect on (7/4 comes from the three-event correlation, we obtain 



14 



and this result is verified with the simulation data as shown in Fig. [TOl 

The covariance between J2 and J3 shows a plateau region for the small value of 
e < 1 and this occurance is due to the fact that J2 and J3 are causally con^elated and 
also that they share a common source of noise. {J2J3) are estimated by considering 
two cases of event sequences: one event of vi comes first and then V2 later, and 
these events occur in the opposite order. By taking into account both the cases, we 
can estimate (J2J3) as 

{■hJs) = ^(^^2(si)i'3(si - 1)> + ^(^^3(sl)^^2(sl)), 
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Figure 10: Flux covariance of two reactions vi and in the two step cascade re- 
action system Fig. 1 A. Parameters: (Xi, X2,pi,p2,P3,P4) = (1,1,0.1,0.01,1,1). 

where the first term represents the case that an event of reaction V2 occurs first, 
resulting in the decrease in Si by one, and then an event of reaction occurs. 
The second term is for the other case that the reactions occur in the opposite order. 
Therefore, we obtain the flux covariance: 

The first term on the left hand side is due to the common source of noise, in this 
case 5*1, and the second due to the causal correlation. The above expression can be 
further simplified to 0-23 ^ ^P2P3{Si). The height of the plateau is well estimated 
at 0.05 (graph is not shown). 

The covariance between J2 and J4 also shows a plateau region for the small 
value of e < 1, and the height of the plateau can be estimated by 

0"24 = {v2iSl)v4,{S2)) - {V2){vi) = P2?'40"i2- 

This estimates the plateau height well (graph is not shown). The reason for the 
occurance of the plateau region is that J2 and J4 have a common source of noise, 
resulting in the flux covariance: E.g., an event of reaction V2 can be correlated with 
that of reaction by events of reaction vi that has occun^ed previously. 

0-34 can be estimated by foUowing the simhar estimation procedure to the one 
for (7^3: 

0-34 =PWi[<yi2 + 2^^-^)] 
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Figure 1 1 : Flux covariance of two reactions vs and in the two step cascade re- 
action system Fig. 1 A. Parameters: (Xi, X2,pi,p2,P3,P4) = (1,1,0.1,0.01,1,1). 

The first term is due to the noise propagation from common sources of noise and 
the second due to the causal correlation. The flux covariance is estimated at 15 (see 
FigHB. 

Finally, for the intermediate and large value of e, i.e., e > 50, four different 
covariance quantities match with one another: cj/3, (T23, (724 , Ji ~ J2 and 
J3 ~ J4. 

In summary, the sum value of the flux CV summation theorem depends on 
which reaction pairs to choose as well as the value of e. The asymptotic forms of 
flux CCVs in the limit of e — > are independent of e, i.e., plateau regions appear, 
if (1) the two reaction steps are affected by the noise propagated from common 
sources or (2) they are directly connected such that one reaction event leads to the 
direct change in the probability that the other reaction occurs. 

3 Estimation of control coefficients for concentration CVs 
from the Lyapunov equation 

In this section, we will show how we estimate control coefficients for concentration 
CV based on the linear noise approximation. Let us define a mathematical notation: 
The matrix component of is ^ denotes the change in the system 
variable x from one stationary state to another due to a parameter perturbation of 

Pi^ Pi + Spi. 

Consider an infinitesimal perturbation in the control parameters denoted by p. 
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The Lyapunov equation Eq.[7] (14 ll. 1421) is invariant because we consider stationary 
state perturbations: 

6{Ja + cr^ + D)=0. 

We obtain 

J— + —J^ + —a + a^ + — = Q, (14) 

6pi 6pi 6pi 6pi 6pi 

where we have used a = <t^ and 6a / 6pi means the change in the concentration 
covariance matrix due to the change in pi and it is defined as an unsealed control 
coefficient of the concentration covariance matrix, a. This unsealed control coeffi- 
cient will be estimated first and then the scaled control coefficient of concentration 
CV will be obtained later. 

To solve the above equation for 6a / 6pi, we need to express 6 J / 6pi and 6D / 6pi 
in terms of concentrations s and p. 6J{s,p)/6pi can be expressed as follows: 

6J{s,p) dJ dJ ds 
6pi dpi ds dpi ' 

By performing the similar procedure for D, 6D/6pi can be expressed as: 

6D{s,p) _ dD dD ds 
6pi dpi ds dpi' 

By substituting the above two expressions in Eq.(fT4l). the unsealed control coef- 
ficients for a concentration covariance matrix (6a/6pi) can be numerically esti- 
mated. 

Next, we need to obtain the control coefficients for concentration CV/CCV 
instead of concentration variance/covariance. The concentration CV is defined as 
Vj^f^ = Ujk/sjSk- The unsealed control coefficients for the concentration CV can 
be obtained: 

6Vjk _ 1 6ajk (Tjk 6sj ajk 6sk 



6pi SjSk 6pi s'jSk 6pi Sjsl 6pi ' 



(15) 



where 6sj/6pi is an unsealed control coefficient for mean concentration sj. In this 
section, we have obtained the mathematical forms of control coefficients, under 
the assumption of the linear noise approximation. At this approximation level, 
the mean concentration dynamics are described by the deterministic rate laws, by 
neglecting the contributions of all concentration covariances and higher moments 



(|43h . Thus, we can express the unsealed concentration control coefficients as in the 
deterministic case: 

6s _ 1 dv 
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where is a reduced stoichiometry matrix (I44t) . By substituting both this ex- 
pression for ds/6pi and the numerical estimate of 6a/6pi in Eq.(fT5]). the unsealed 
control coefficients for the concentration CV/CCV can be estimated. 

Finally, we convert the unsealed control coefficient to a scaled version by using: 



Cn 



We provide a MATHEMATIC A file for this estimation in the Supporting MATH- 
EMATICA file. 
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